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Abstract. Einstein's field equations for timelike self-similar spherically 
symmetric perfect-fluid models are investigated. The field equations are rewritten 
as a first-order system of autonomous differential equations. Dimensionless 
variables are chosen in such a way that the number of equations in the coupled 
system is reduced as far as possible and so that the reduced phase space becomes 
compact and regular. The system is subsequently analysed qualitatively using the 
theory of dynamical systems. 

Using this approach, we obtain a clear picture of the full phase space and the 
full space of solutions. Solutions of physical interest, e.g. the solution associated 
with criticality in black hole formation, are easily singled out. We also discuss 
the various 'band structures' that are associated with certain one-parameter sets 
of solutions. 
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1. Introduction 

Spherically symmetric perfect-fluid models admitting a homothetic Killing vector has 
attracted considerable attention during the last couple of decades (see, e.g., |Q, |j and 
references therein). The main reason for this is that the homothetic Killing vector 
reduces the field equations of spherically symmetric perfect-fluid models from partial to 
ordinary differential equations. This simplification does not mean that these so-called 
self-similar models do not contain interesting physical phenomena - quite the contrary. 
The possibility of violation of the cosmic censorship hypothesis Q and the existence of 
sonic surfaces (shock waves) D are some examples. The growth of primordial black 
holes and the evolution of voids can also be studied using self-similar perturbations 
of the flat Friedmann-Lemaitre-Robertson- Walker (FLRW) model ^, 0] . In recent 
years it has become clear that self-similar models are of crucial importance for the 
understanding of black hole formation through gravitational collapse. In the full 
space of spherically symmetric solutions, certain self-similar solutions are separatrices 
between dispersive solutions and solutions developing black holes p|, |lOt. 

Self-similar spherically symmetric perfect-fluid models exhibit several preferred 
geometric structures. There are two preferred directions corresponding to the 4- 
velocity of the fluid and the homothetic Killing vector, respectively. Adapting the 
coordinates to the 4-velocity of the fluid leads to the 'comoving fluid approach' 
[ PI H, ^ |l2|, Adapting the coordinates to the homothetic Killing vector leads 
to the 'homothetic approach' |l4|. The spherical symmetry surfaces also constitute 
a preferred geometric structure. The area of the spherical symmetry surfaces is used 
as a coordinate in the 'Schwarzschild approach' (see, e.g., |3|). There are also some 
other possibilities, e.g. synchronous coordinate systems and null coordinates, which 
are useful for discussing singularities All these approaches are complementary 
and a full understanding probably relies on a combination of results obtained in the 
various approaches. 

In this paper we will use the diagonal homothetic approach. This has the 
advantage of yielding equations which are very similar to those of hypersurface 
homogeneous models. There exists a wealth of literature on how to treat such models. 
Thus many of the ideas obtained in this area can be carried over into the realm of 
self-similar spherically symmetric models. However, there are of course disadvantages 
as well. The symmetry surfaces will in general change causality. Thus, in the diagonal 
homothetic approach the spacetime has to be covered with two coordinate patches; 
one for when the homothetic Killing vector is spacelike and one for when it is timelike. 
Then these two regions must be joined where the homothetic Killing vector is null. 

In a previous paper IT^ we considered the spatially self-similar (SSS) spherically 
symmetric models. In the present paper we consider the case of timelike symmetry 
surfaces, constituting the so-called timelike self-similar (TSS) spherically symmetric 
models. We obtain a full phase space picture which leads to a more complete 
understanding of these models than has been previously obtained. Thus the 
work presented here confirms, corrects and expands the analysis of several authors 

The TSS spherically symmetric models are characterized by a four-dimensional 
homothetic symmetry group acting multiply transitively on three-dimensional 
timelike hypersurfaces. The line element, written in diagonal form where one of the 
coordinates is adapted to the homothetic symmetry, takes the form [Q 

ds^ = e^'ds^ = e^* [-Dl{x)dt^ + dx^ + Dl{x) {d9^ + sin^ edip^)] . (1) 



Timelike self-similar spherically symmetric perfect-fluid models 



3 



We will consider perfect-fluid models. The energy momentum tensor, Tab, is given by 

fab = fiUaUb + p{UaUb + Qab), (2) 

where fi is the energy density, p is the pressure, and the 4-velocity of the fluid. We 
will assume 

P=(7-1)A, (3) 
as an equation of state where the parameter 7 takes values in the interval 1 < 7 < 2, 
which include radiation (7 = |). We have excluded dust (7 = 1) and stiff fluids 
(7 = 2) as they behave quite differently compared to those in the above interval and 
need special treatment. The dust solutions are explicitly known (these models are just 
special cases of the general spherically symmetric dust solutions, which arc all known 
explicitly, see, e.g., 0, 

The outline of the paper is the following. In section 2 Einstein's field equations 
are rewritten in terms of a dimensionless set of variables, in order to obtain a maximal 
reduction of the coupled system of ordinary differential equations. The variables are 
chosen so that they take values in a compact phase space. In section 3 the equations 
are subsequently analysed by using the theory of dynamical systems. In section 4 a 
numerical investigation is undertaken and global dynamical features are considered. 
Appendix A describes the properties of the fluid congruence and also the condition 
for the spacetimes to belong to Petrov type 0. Appendix B gives the relation between 
various coordinates and the relation between the present variables and other variables 
used in the literature. 



2. The dynamical system 

Following we introduce 

D, = Br' - e^''-'^^ D, = B,-' = e^'^^^ (4) 

e =3/3", cr+=3^+, 
where a dot denotes d/dx. The quantities 9 and a+ describe the kinematical properties 
of the normal congruence of the symmetry surfaces in the static (M, ds^) spacetime 
that is conformally related to the physical spacetime (M, ds^) with the homothetic 
factor e^* (see, e.g., jl^); is the expansion while a+ describes the shear of the 
normal congruence of the symmetry surfaces in (M, ds^). Expressed in an orthonormal 
frame, the fluid 4-velocity is parametrized by = (1, u, 0, 0) / \/l — u^. Einstein's field 
equations, Gab — Tab, and the equations of motion for the fiuid, T'°^;fc = 0, lead to a 
set of evolution equations for 9, a+, Bi, B2 and u. In addition, a constraint equation, 
and equations for /it and /it are obtained. The quantity /tt, related to fi by 



1 + (7~ 1) ^^ 

1 -u2 



= \ 5 e jl, (5) 



is the energy density of the fluid as measured by an observer with a 4-velocity 
associated with the homothetic symmetry. 

As in the SSS case we will now choose to use the 0, (t+ variables, defined by 

e^^{2e-cy+), a+ = ^{-d + 2a+), (6) 

This is done in order to simplify the constraint, while still keeping the quadratic form 
of the defining equation for /tj (equation (g) below). This leads to: 
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£7+ = 

Bi = 
B2 = 
u = 

Constraint equation 



1 

1 

75 



0-+ + 9a+ - 3Bt - 3^ 



(7 -l)(l-^^) 
+ (7 - 1)^2 



+ 29a+ + 



3 (37 - 2) + (2 - j)i 
2 1 + (7-1)^2 



-Ait 



1 



(7) 



V37K-(,-l)]{-[^(^-^)^" + ^^-] 
+%/3 [(7 - 1)(37 - 2) - (2 - 7)«2] Bi} 



7«Mt - [1 + (7 - ^+Bi = 0. 



Defining equation for jit 



Mt = 



1 + (7 - 1)m2 
m2 + (7 - 1) 



[^2 - - Wl - 3S|] 



(8) 



(9) 



Auxiliary equation 



3 [1 + (7 -1)^*2] [u2_(^_l)] 



pj.{V3[(40>a+)7-60>^ 



+ 6(-77 + 372 + 4)wBi 
-73(7- l)(2^ + a+) 



4 

u 7 ■ 



l(-57 + 272 +4)u3Bi 
a/3ct+7| . 



(10) 



The defining equation for /it shows, assuming non-negative energy density, that 9 is 
a 'dominant quantity'. Consequently, in order to obtain bounded 'reduced' variables, 
we now introduce ^-normalized dimensionless variables S+, A, K 



<7+ 



A: 



VSBi 



e 



K = 3 



B2 
9 



The density fit is replaced by the density parameter Qt, which is defined by 

* 6(2 

The introduction of a dimensionless independent variable t] 

dx _ 

d^~~f' 
leads to a decoupling of the ^-equation 

(7-l)(l-«2) 



(11) 
(12) 

(13) 



= -{l + q)9, 



E+(l + E+)-A^ 



(14) 



l + (7-l)u2 

where a prime denotes d/drj. The remaining coupled evolution equations can now be 
written in dimensionless form: 
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Evolution equations 

Si : 



1 - si + 



1(37 -2) + (2 -7)1*2 



(7-l)(l-^^) 
1 + (7-1)1*2 



2 1 + (7 -1)1*2 



A' = 
K' ^2 



1 



2S+ + - - ^ 



(7-l)(l-"') 



+ (7 - 1)1*2 

2> 



A, 



(15) 



l + (7-l)u2 



if, 



1-U2 



7[m2 - (7- 1)] 

{7 [2(7 - 1) + 7S+] u + [(7 - 1)(37 - 2) - (2 - 7)^2] A] , 



Constraint equation 

G = "fuQt - 2 (1 + (7 - 1)m2) t+A = 0. 



Defining equation for fit 

^ 1 + (7 - 1)^^" 

* 1/2 + (7 - 1) 



(16) 



(17) 



Auxiliary equation 



^t 



S+7 



(18) 



[1 + (7 - l)l/2] [l/2 _ (-^ _ 1)] 

{ [(4 + i;+)7 - 6] u^j + 2(-77 + 87^ + 4:)uA 
- 2(-57 + 272 + 4)u^A - (7 - 1)(2 + S+)?/7 - 
+2(1 + g) [1 + (7 -1)1.2] [^2_(^_^)]|_ 

Equation (||) makes it impossible for 9 to change sign when /if 7^ 0. If ^ > 0, then 
A > 0, where we have included the zero boundary value. The field equations are 
invariant under the transformation 

it+,A,K,u)^it+,-A,K,-u). (19) 

Thus it is trivial to obtain the case when 9 < corresponding to A < 0, once the case 
A > has been analysed. We will therefore assume that ^ > 0, A > 0. Note that 
both cases are needed in order to obtain a global picture. 

The surfaces defined by li^ = 7 — 1 are surfaces of non-extendibility of solutions. 
The only way to cross such a sonic surface analytically is through the line where 
the numerator of the u' equation of ( p^ vanishes. This is the sonic line. To study 
the behaviour of trajectories around this line, we make the following non-monotonic 
change to a new independent variable H, 

d?7 ~ (7- 1) - 1/2- ^ ' 

The line element can be obtained when K, A and 9 have been found through the 
relations 



Di = V^i9A)-\ D2 = Vi{9^K)-^ 



/2 



a/3 



drj 



(21) 
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Table 1. 


Equilibrium points of the TSS phase 


space. 
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Notation 
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o 7-1 


n 7^+4(7-1) 
" (37-2)^ 





4(7-1) 
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M± 





1 


see text 





SL 




see text 


-V7- 1 






1 





±1 







-1 





±1 





M+ 





1 


1 
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1 


-1 
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S+ + 1 


-1 









1 


±1 






Table 2. Linear analysis of equilibrium points of the TSS phase space. 

Eigenvalues 



Notation 


Elim. 


Ai 


A2 


A3 A4 


K 

CO 


A 
A 


V7-6 
7-1 
2-7 

7-1 

1 


2 
2 
2 


37-2 

7-1 
2-7 
7-1 

-1 -2 


T 

M± 


A 
K 


2-7 
37-2 

-2 


see text 


see text 


SL 


K 











A 
A 


4 



2 
2 


2(37-2) 
2-7 

2 


M+ 
M- 


A 
A 


-4 



-2 

-2 


2(57-6) 
2-7 

-2 


H- 


A 





-2(1 + 22+) 


-2(1 + 2S+) 


C± 


K 


-1 


1 


4(7-1) 
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3. Dynamical systems analysis of the reduced phase space 

The reduced phase space is determined by X,?/), related by the constraint 

G = 0, given by equation (p^). As in JTst , we include the boundary in order to 
obtain a compact phase space. The boundary is given by a number of invariant 
submanifolds: Vtt = Q, A = {), K = and m = ±1. The constraint cannot be solved 
globally everywhere. Instead we will follow ^ and solve it locally around the 
equilibrium points (i.e. we will solve the linearized constraint for different variables 
at different equilibrium points). This formulation will enable us to achieve a good 
understanding of the global structure of the dynamics of the reduced phase space. For 
an introduction to dynamical systems analysis see, e.g., [^, ch 4. 

3. 1 . Equilibrium points and local analysis 

There are numerous equilibrium points of the TSS dynamical system. In table 1 the 
equilibrium points are presented together with f2t, which will indicate if a point is 
on the vacuum boundary fij = or not. Around each equilibrium point we locally 
eliminate one of the variables by solving the constraint to linear order. In table 2, 
the variables eliminated in the local analysis are listed together with the eigenvalues 
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at each point. The equilibrium points are often closely related. Throughout, they 
are denoted as Kernef ^"fe"' ^ . When there is no risk for confusion we have suppressed 

sgn(S+) or sgn(u). The kernel indicates the interpretation of the point: M and C 
refer to Minkowski spacetime; K indicates a Kasner solution; F is the flat FLRW 
solution and T refers to a static solution discussed below. Note that the self-similar 
solutions under investigation correspond to orbits in the interior of phase space that 
asymptotically approach the various equilibrium points. Below, we will comment on 
some of the equilibrium points. 



The equilibrium point C°. This point corresponds to the Minkowski spacetime 
expressed in spherically symmetric coordinates. The constraint surface is degenerate, 
i.e. VG — (0,0,0,0). Hence we keep all four eigenvalues. 



The equilibrium point T . The eigenvalues A2 and A3 are 



^^^^-"2^* 2(3^^2^ ■ 

Note that these eigenvalues are always complex for 1 < 7 < 2, and that the real part 
is always negative. The point T corresponds to a self-similar static solution associated 
with many names: Tolman Oppenheimer and Volkoff |Q, Klein and Misner 
and Zapolsky |Q. The same solution also leads to an orbit in the interior phase space, 
corresponding to a different foliation of the spacetime. This orbit will be referred to 



as the static orbitnl, and will be discussed in section 4.2 below 



The equilibrium points . For these points, the variable u takes the values 



7(7 - 1) ± - 1) (7^(7 - 1) + (37 - 2)(2 - 7)) 

u= ^ . (23) 

2-7 

The expressions for A2 and A3 will not be given since they are rather complicated. 
Instead we briefly comment on the stability of the points. At 7 = 1 both points 
coincide, but as 7 increases the points move apart. For 7 = | the point Af+ passes 
through the point Af + and leaves the physical part of the TSS phase space. It should 
be noted that the point A/+ leaves the TSS phase space only to appear in the SSS 
phase space as one of the points ±M The point M~ exists for all values of 7 in 
the interval 1 < 7 < 2. It is important to note that M~ always has — -^7 — 1 < u < 0, 
while M"*" always has u > -^7 — 1. This behaviour can be seen in figure 2 (c) - (f). All 
eigenvalues are negative in the interval 1 < 7 < | for M"*". Thus this point is always 
an attractive node. For two eigenvalues are negative and the third is positive in 
the interval 1 < 7 < 2. Hence this point is a saddle with one outgoing eigendirection. 



The equilibrium points SL. The flows are opposite to each other on the different sides 
of the sonic surfaces at = 7 — 1. Thus it follows that orbits, in general, cannot pass 
through these surfaces continuously. It is only possible to cross a sonic surface along 
the sonic lines S'L, where the numerator of the u' equation in (15) vanishes. These 
lines are defined by 

(,) 7(2(7-l)+7S+) ^_(+)^/— T f24) 
^ - + 4(7-1)3/2 ' M - _ V7 I24j 

X In (li, it was called the 'TOVKMZ' orbit. 
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where the lines of equilibrium points are parametrized by The signs within 

parentheses correspond to the line for which u = — 1. This line of equilibrium 
points is located outside the physical phase space for 1 < 7 < 2. Consequently the only 
relevant sonic line is the one at u = —\/"f — 1- Using the 'shock-adapted' independent 
variable defined in yields the eigenvalues A2 and A3, which are solutions to the 
equation 



where 



4(7-l)A' + (2-7)(6A + c) = 0, (25) 

b = -27 (2(7-1) +7S+), c - (2 - 7) [ciE^ + C2S+ + C3] , 

01 = 77^-3872 + 867-8, C2 = 8(7- 1)(272 - 77 + 4), (26) 

C3= -4(2-7)(7-l)2. 

An analysis of the stability of the line shows that it splits into four parts, see figure 1. 
The boundaries of the physical interval of the line are determined by A ^ and 
K = 0. The line is a saddle at ^ = 0, then it turns into an attractive node 



containing the static orbit, discussed in section 4.2. This attractive nodal region 
is called the S.a.n. region. The eigendirection associated to the eigenvalue with 
the smaller absolute value dominates the dynamics close to the sonic line 10. This 
direction is called the primary eigendirection. The other direction is referred to as 
the secondary eigendirection. Subsequently, the eigenvalues become complex and the 
stability changes to an attractive focus. Orbits approaching the sonic line in this region 
will spiral and hit the sonic surface infinitely many times, and are thus unphysical 
Finally, the stability changes back to an attractive node, now containing the flat 



FLRW orbit, discussed in section 4.2. This attractive nodal region will be called the 
F.a.n. region. The behaviour along the sonic line is similar for all values of 7 in the 
interval 1 < 7 < 2. Note that the stability is with respect to the new independent 
variable defined in (^). In the following, if not explicitly pointed out, all figures and 
discussions of phase space fiows are with respect to the original independent variable. 
Note that there appears to be an error in the stability analysis of Bogoyavlensky, 
which leads to an erroneous picture of the flow on the phase space boundary Q . At 
the boundaries of the focal region, the two non-zero eigenvalues are equal. These 
degenerate nodes will be denoted S.d.n. and F.d.n., respectively. For 7 = |, both 
the static orbit and the flat FLRW orbit coincide with the corresponding degenerate 
nodes on cither side of the focal region. 



The equilibrium point AI^ . Note the change in stability of this point at 7 = |. This 
bifurcation corresponds to when the point A/+ passes through M+, leaves the TSS 
phase space and enters the SSS phase space. 



The equilibrium points TL^ . This line of equilibrium points is an artifact of the 
diagonal homothetic approach. It is associated with a causality change of the 
homothetic vector field. To obtain a global picture, one needs to match the TSS phase 
space with the SSS phase space. For the the choice ^ > (0 < 0), corresponds to 
the line of equilibrium points +7i^ (-7i^) in the SSS phase space JlSt . 
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Figure 1. This figure shows the stability of the sonic fine, parametrized by S-|_, 
for 7 in the interval 1 < 7 < 2. The attractive nodal regions are denoted by F.a.n. 
and S.a.n., while the attractive focus is denoted by a.f. The full curves mark the 
boundaries of the physical interval for S4.. The dotted curves are the values on 
the sonic line for the flat FLRW orbit and the static orbit, respectively. A similar 



3.2. Invariant suhmanifolds on the boundary of the TSS phase space 

As stated previously, the boundary is described by a number of invariant subsets: 
= 0,A = 0,K = 0,u = ±1. The constraint @ leads to uQt = when ^ = 0. 
For u = one obtains the reduced equations for the static spherically symmetric 
perfect-fluid models. When fit = 0, u 7^ one obtains the static vacuum equations for 
a spherically symmetric space time with a test fluid. The submanifolds u — ±1 are 
described by the same equations as the TSS spherically symmetric model with a null 
(e.g. neutrino) fluid. The equations for the submanifold ^ = are the same as the 
reduced equations for the type il models. The submanifolds are summarized in table 
3 where we also introduce designations. 

We will now describe the dynamical features of the individual boundary 
submanifolds. 

The submanifolds. These submanifolds correspond to Cauchy horizons that 
are black hole event horizons or, in a cosmological context, particle horizons Q. 




Table 3. The various boundary submanifolds. 
Boundary Restriction 



u = ±l 
V S-|-=0, n = o 

il ^" = 

SV± A = 0,n = 0, sgn(M) = ±1 

S A = 0,u = 
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Figure 2. The phase portraits of the invariant submanifolds that constitute the 
boundary of the reduced phase space. 

The dynamical structure of the N~ submanifold is given in figure 2(a) and the 
corresponding diagram for 7V+ is given in figure 2(b). The 7V^ submanifolds are 
solvable. The integral describing the various orbits is given by 

z ^-o = constant. (27) 

[l-S+-sgn(«)A]' 

The vacuum V submanifold. The dynamical structure of the V submanifold is given 
in figure 2(c) (7 < |) and 2(d) (7 > |). Note that A is an increasing monotonic 
function, which by the constraint implies that ^ is a decreasing monotonic function. 
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The il submanifold. The dynamical structure of the plane-symmetric type 
il submanifold is given in figure 2(e) (7 < |) and 2(f) (7 > |). 

The static vacuum SV± submanifolds. The dynamical structure of the SV^ subman- 
ifold is given in figure 2(g) and the corresponding diagram for SV+ is given in figure 
2(i). The SV± submanifolds are solvable. The integral describing the various orbits is 
given by 



(1 + S+)2-7E4(^- 



= constant. (28) 



The static S submanifold. The dynamical structure of the S submanifold is given in 
figure 2(h). There exists a monotonic function for this submanifold, excluding the 
point T. It can be found by using the Hamiltonian methods developed in ch 10 of 
[^ . It is given by 

^ (37- 2 + 2(7- 1)S+)' ^ 4(7-1)^ ^ (7^+4(7-1)) 

npRi (57-4)7''^ (57-4)7 ' 



with 



Z' _ 2(37-2)[2(7-l)-f-(37-2)S 



1 2 



+ 



Z [37 - 2 + 2(7 - 1)S+J (57 - 4)7 ■ ^^^^ 

This monotonic function prevents the existence of equilibrium points, periodic orbits, 
recurrent orbits and homoclinic orbits in this region, see e.g. The static 

spherically symmetric perfcct-fiuid models have been studied qualitatively by Collins 
using non-compact variables but the above monotonic function is, to our 

knowledge, new. The compact formulation taken together with this monotonic 
function leads to a complete understanding of the dynamics of static spherically 
symmetric models. It is worth noting that the orbit from C° to T corresponds to 
the only regular solution which exists for these models (modulo a scale invariance 



parameter). As will be discussed in section 4.2, this orbit constitutes part of the 
boundary of a two-dimensional interior submanifold corresponding to regular self- 
similar solutions. The mass function (see Appendix A.2| ) for the S submanifold is 



^^l-ii±3:l!. (31) 
R K ^ ' 

The condition m = defines a curve between and C". This curve is indicated by 
the dotted curve in figure 2(h). Both the equilibrium point T and the regular solution 
C^-T are located in the m > region. The corresponding solutions are the only static 
solutions which have a non-negative mass everywhere, since all other orbits come from 
K^, located in the m < region. 

4. Global behaviour 

By appropriately 'gluing together' the boundary submanifolds of the previous section, 
we obtain the reduced phase space shown in figures 3 and 4. The TSS perfect-fluid 
models correspond to orbits in the interior reduced phase space. The il submanifold 
makes up the bottom of a 'tent' with the other submanifolds as 'walls'. The point C'^ is 
the top of the 'tent'. Indicating the stability of the equilibrium points on the boundary, 
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Figure 3. This picture sliows how the boundary of the reduced phase space of 
the timelike self-similar spherically symmetric models is constructed out of the 
invariant submanifolds discussed in section |3.2| . The broken curves describe the 
sonic surfaces. 



we obtain figure 4 for 7 > |. For 7 < |, the point Af+ enters the physical phase 
space, but is always beyond the sonic surface located at u = -^7 — 1. It is impossible 
to analytically continue solutions through the sonic surface there, in contrast to the 
one at u = — -^7 — 1 . Thus the bifurcation at 7 = | only affects the 'physically 
uninteresting' part of the phase space, ^ u < 1. The advantage of a compact 

and regular phase space is now apparent: no parts of the phase space are 'crushed'. 
This is in contrast to, e.g., Bogoyavlensky Q, where parts of the phase space are cut 
off, while others are located at infinity. 

The orbits on the boundary corresponding to the eigenvector directions of the 
points M~ and C° are shown in figure 4, as is the static orbit along the eigenvector 
direction of the point T pointing into the phase space. The two-dimensional separatrix 
surface entering the interior phase space from C° and bounded by these orbits will be 



discussed in section 4.2 



4.I. Monotonic functions 

As indicated above, monotonic functions are important tools for understanding the 



global dynamics. The function |14 



whose derivative with respect to the independent variable 77 is 
(2-7)(37-2)(l-u2)A] 



2\-(2-7) 



F' 



F, 



(32) 



(33) 



is easily seen to be monotonic in the regions u < and m > 0. This monotonic function 
tells us that all orbits in the m > part of the phase space come from the sonic surface 
at u = -^7 — 1 H|. Thus, all the 'interesting' dynamics takes place for u <0. 
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Figure 4. The global structure of the reduced phase space for 7 > ^. The orbits 
on the boundary corresponding to the eigenvector directions of the points M~ 
and C*^ are shown, as is the static orbit along the eigenvector direction of the 
point T pointing into the phase space. 



4.2. Invariant suhmanifolds in the interior of the TSS phase space 

A number of trajectories, corresponding to exact solutions of the field equations, can 
be found. 

The flat FLRW orbit. An orbit corresponding to the flat FLRW solution can be 



obtained by imposing the Pctrov type condition, C = see equation (A.4), along 
with vanishing fluid shear. It is described by 

E A- f34) 

+ (37-2) + 4u2' (37-2) + 4u2- ^ ' 

The flat FLRW orbit passes through the sonic line at an equilibrium point in the 
F.a.n. nodal region. The value of the variable E+ on the sonic line is given by the 
upper dotted curve in figure 1. The flat FLRW orbit approaches the sonic line along 
an eigendirection. Thus the flat FLRW solution is analytic across the sonic surface, 
as expected (for a discussion about differentiability, see the discussion about the ASL 
submanifold below). The corresponding eigendirection is primary for 7 < |, and 
secondary for 7 > | ■ This was pointed out in |l^ . 

The static orbit. The static orbit, discussed in a previous section, enters the physical 
phase space from the point T and can be found explicitly: 

- _ 2(7-1) ^_ -ju 

Note that the static solution is not the only solution that is represented both as an 
equilibrium point on the boundary and as a trajectory in the interior phase space. This 
is also the case with the flat FLRW solution in the SSS sector, see [^. The static 
orbit passes through the sonic line at an equilibrium point in the S.a.n. nodal region. 
The value of the variable E+ on the sonic line is given by the lower dotted curve in 
figure 1. The static orbit approaches the sonic line along an eigendirection. Thus the 
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static solution is analytic across the sonic surface, as expected. The corresponding 
eigendirection is secondary for 7 < |, and primary for 7 > |. 

The RC submanifold. A requirement on solutions describing collapse is that they 
should initially be regular at the origin. The point C*^ is associated with such a 
regular centre ^, |l^. The outgoing eigendirections at this point span a surface 
which enters the reduced phase space. This one-parameter family of orbits will be 
denoted the regular centre or RC submanifold. Its boundary on the V submanifold is 
the orbit from C", terminating at M~ (see figure 2(c) and (d)). The other boundary 
is the orbit from to T, continuing into the phase space along the static orbit. The 
orbits in the RC submanifold can be parametrized by the density at the centre. For 
example, Ori and Piran use a parameter Dq (see Appendix B.3). 

The closer an orbit starts to the — T orbit the more it spirals around the static 
orbit. This is to be expected, since the T equilibrium point has a pair of complex 
eigenvalues. The circulation around the static orbit is associated with a sign change 



of the radial 3- velocity of matter vji (see the expression in Appendix B.4). For orbits 
close to the static orbit, vr changes sign several times. This corresponds to several 
collapsing and expanding spacetime regions enclosing each other. Thus solutions fall 
into different physical families depending on the number of zeroes of vji. This 
defines a band structure of the RC submanifold. Near the V submanifold we have 
'pure collapse' solutions (n„ = 0)- It is obvious that there is an infinite number 
of n„-bands (since T has a pair of complex eigenvalues). The flat FLRW solution 
corresponds to an orbit within the RC submanifold, and belongs to the ti^, = riy- 
band. The higher n„-band structure corresponds to orbits close to the — T orbit. 
These orbits and their continuation along the static orbit correspond to spacetimes 
that can be described in terms of self-similar coordinates as follows: the spacetime 
(the subsonic part) can be divided into two regions. In the first region the solution 
is approximated by the regular static solution multiplied by the homothetic factor. 
Thus the solution is approximately conformally static in this region. In the second 
region the solution is approximated by a perturbation of the static self-similar solution 
(written in conformally static form). This second regime thus corresponds to 'quasi- 
static' oscillations around the static solution. The static orbit itself is the only orbit 
starting at T, which implies that there are no asymptotically quasi-static solutions. 
However, there is a two-parameter set of solutions coming arbitrarily close to the static 
solution during a limited part of their evolution. 

The intersection of the RC submanifold with the sonic surface exhibits another 
band structure in the following sense: when varying the density parameter it is 
found that there are continuous intervals of Dq for which the orbits hit the sonic surface 
at the sonic line. The stability of the sonic line determines if an orbit correspond to 
a physically interesting solution. Consequently, 'physical' orbits must hit the sonic 
line in one of the nodal regions, or approach the sonic line along the attractive 
eigendirection of the saddle region. These bands were considered by Ori and Piran 
[^, and we will denote them OP-bands. One OP-band consists of orbits hitting 
the F.a.n. nodal region of the sonic line. This is the first band in Ori and Piran's 
terminology. For 1 < 7 ^ 1.04, its upper boundary is an orbit approaching the sonic 
line along a secondary eigendirection of the F.a.n. nodal region. The corresponding 
solution is known as the general relativistic Penston-Larson (GRPL) solution For 
1-04 ^ 7 < |, the boundary is the F.d.n. degenerate node. This is in accord with 
the proposition that orbits in the neighbourhood of the flat FLRW orbit meet the 
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log(Do) 



Figure 5. Values of S_|_ at the sonic surface for orbits with 7 = 1.4. The broken 
line is the static value, and the dotted line is the S.d.n. degenerate node. One OP- 
band starts near the saddle region for log(Do) ~ 0.81 (not shown) and ends when 
the curve moves into the focal region at log(Do) ~ 1.05. This OP-band is the 
second band, using Ori and Piran's terminology. The third OP-band begins when 
the curve moves into the S.a.n. nodal region at log(Do) ~ 1.5. This OP-band is 
infinitely broad, as the curve never again moves into the focal region. 



sonic surface along the sonic line for 7 < | 14 1. For higher 7's, the boundary is the 
flat FLRW orbit. The first OP-band is separated from higher OP-bands by a broad 
region where the orbits hit the sonic surface at other places than the sonic fine. Orbits 
spiraling around the static orbit oscillate around the static value when they hit the 
sonic line. As discussed in connection with figure 1, the static orbit coincides with 
one of the degenerate nodes for 7 = | • This implies that there will be a large number 
of OP-bands for equations of state with 7 near |, because each time solutions move 
into the focal region on the sonic line, an OP-band ends with the S.d.n. degenerate 
node as boundary. When solutions move back into the S.a.n. nodal region, the next 
OP-band begins, again with the S.d.n. degenerate node as boundary, see figure 5. 
As Do increases, the amplitude of the oscillations decreases. Consequently, there is a 
value Dod such that Do > Dq^ implies that orbits remain in the static nodal region. 
Thus, there is a finite number of OP-bands, except for 7 = | for which the static 
orbit coincides with the S.d.n. degenerate node. Note that the OP-band structure 
also is affected by whether the orbits hit the sonic surface on the sonic line or not. 
This is the dominant effect for creating OP-band structure for soft equations of state. 
For 1 < 7 < 1.41, the lower boundary of the second OP-band corresponds to an 
attractive eigendirection in the saddle region. For 1.41 < 7 < 1.89 the boundary is a 
secondary eigendirection in the S.a.n. nodal region, while for 7 > 1.89 it is the S.d.n. 
degenerate node (see also the discussion about critical behaviour below). The upper 
boundary changes character, depending on the value of 7. For 1<7<1.11, itisa 
secondary eigendirection in the S.a.n. nodal region. For 1.11 ^ 7 ^ 1.45, it is the 
S.d.n. degenerate node. The orbits in the second OP-band belong to n,u-band = 1 



Timelike self-similar spherically symmetric perfect-fluid models 



16 



or — 2. For 7 > 1.45, the orbits never move into the focal region on the sonic 
line. The OP-band structure then degenerates to only two bands; the first OP-band, 
containing the flat FLRW solution, and a second band, infinitely broad in Dq. 

The ASL submanifold. In order for the solution to be analytic (or C°°) at the sonic 
point, the corresponding orbit has to approach the sonic line along an eigendirection 
in one of the nodal regions ||^, or possibly along the attractive eigendirection in 
the saddle region. The corresponding orbits form a submanifold which we will call the 
ASL submanifold. 

The RCASL submanifold. An analytic solution is regular at the centre and analytic 
at the sonic point. These solutions correspond to orbits that belong to the intersection 
of the RC and ASL submanifolds. This intersection forms a submanifold which we 
will denote as the RCASL submanifold. The RCASL submanifold turns out to be a 
discrete set of orbits. Typically, there is only one such solution for each n„-band. A 
convenient numerical way to find orbits belonging to RCASL — RC n ASL is to start 
along eigendirections at the sonic line. This gives a set of orbits in ASL. For an orbit 
belonging to RC, should tend to zero, as = at C". It is thus possible to 'fork 
in' the interesting solutions. 

Foglizzo and Henriksen [ p^ use a set of variables particularly adapted to numerical 
studies of orbits between the sonic surface and the regular centre. However, we point 
out that their variables are non-compact and distort the global picture. 

Criticality and self- similarity. Maison ||^ and Hara et al | p^ , p7| consider the critical 
behaviour of spherically symmetric perfect-fluid collapse. A particular self-similar 
solution turns out to be of fundamental interest in these studies. This solution 
corresponds to an analytic orbit in the n„ = 1 n^-band, and constitutes the lower 
boundary with respect to Dq of the second OP-band for 1 < 7 < 1.89. For 7 < 1.41 it 
approaches the sonic line in the saddle region along the attractive eigendirection. 
In the interval 1.41 ^ 7 < 1.89, it approaches the sonic line along a secondary 
eigendirection in the S.a.n. nodal region. For 7 « 1.61, it passes through the point 
for which the static orbit reaches the sonic line. Above this value of 7, the zero in vr 
occurs on the opposite side of the sonic point. For 7 > 1.89, the solution ends in the 
focal region, and thus becomes unphysical, see figure 6. The S.d.n. degenerate node 
is then the lower boundary of the second OP-band. 

The F submanifold. For a global picture of self-similar spherically symmetric perfect 
fluid models, a careful study of the matching between the TSS and the SSS regions 
is necessary. In our previous work on the SSS region of the phase space we 
studied a separatrix surface called the F submanifold. Orbits in this submanifold 
originate from an equilibrium point -^F in the SSS region of the phase space. This 
point corresponds to the flat FLRW solution and it gives the submanifold its name. 
However, the flat FLRW solution is not only represented as an equilibrium point on 
the boundary, but also as an orbit within the F submanifold (corresponding to a 
different foliation of the spacetime). The remaining orbits within the F submanifold 
correspond to models that are of considerable physical interest since they can be 
interpreted as density perturbations of the fiat FLRW model. The F submanifold in 
the SSS region is divided into two parts; one corresponding to orbits lying entirely in 
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Figure 6. This figure shows the stabiUty of the sonic Une, parametrized by 
E+, for 7 in the interval 1 < 7 < 2. The heavy full curve corresponds to the 
intersection of the critical solution with the sonic line. For other designations, see 
figure 1. 



the SSS region and one corresponding to orbits reaching the 7i~ hne of equihbrium 
points (this equihbrium hne is denoted +7i~ in the SSS region). This latter part of 
the F submanifold corresponds to solutions that can be analytically extended into the 
TSS region. 

In order to connect the SSS and TSS phase spaces we use the variables of Foglizzo 
and Henriksen ||l3|| (these variables arise in a comoving context and can therefore be 
used to deal with the situation where the self-similar symmetry surface becomes null, 
see Appendix B). Orbits corresponding to overdense solutions diverge from the flat 
FLRW orbit. Sufficiently overdense solutions re-enter the SSS region through H~ and 
correspond to black hole solutions. However, there is a set of overdense orbits reaching 
the sonic surface. 

Investigating the extendibility of orbits through the sonic surface, we find a band 
structure similar to the OP-band structure of the RC submanifold. We refer to these 
bands as F-bands. All underdense orbits, from the orbit in the il submanifold to the 
orbit corresponding to the flat FLRW solution, pass the sonic surface through the 
sonic line. For 7 > | the limit of this band is the flat FLRW solution. For 7 < | 
there is a set of orbits corresponding to overdense solutions in this band. Between the 
flrst F-band and the next F-band there is a one-parameter subset of orbits that are 
inextendible and end at the sonic surface. This is similar to the region between the 
flrst and second OP-band of the RC submanifold. Between the subset of inextendible 
orbits and the orbits that return to the SSS region, there is a small band of orbits 
corresponding to overdense solutions that extend through the sonic surface via the 
sonic line. Our numerical investigation indicates that there is just one overdense F- 
band, rather than a hierarchy of several overdense bands (as was the case for the 
overdense OP-bands of RC for 7 < 1.45). The interval of the sonic line corresponding 
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to the overdense F-band is wider than the corresponding interval covered by these OP- 
bands. Consequently, there are overdense subsonic {u > — -^7 — 1) solutions which 
are regular at the sonic surface, but cannot be connected with C'^. Numerical analysis 
indicates that these orbits do not recoUapse, but have a second irregular sonic point. 
Therefore, they are unphysical This is in accord with the conjecture of Carr 

and Coley that asymptotically Friedmann solutions which contain black holes are 
supersonic everywhere. 



Mass, energy and density fluctuations. The mass function 2m/ R (see equation (A. 5) 
in Appendix A. 2) is a physically important quantity. By imposing the condition 
m > 0, only certain parts of the phase space are physical. In figure 2(h), the dotted 
line indicates the condition to = in the static S submanifold (see the discussion of 
this submanifold in section 



3.2) 



Another interesting quantity is the asymptotic energy per unit mass £^ of a 
spherical shell (see equation (A. 6) in Appendix A. 2). Carr and Coley ||l|, |^ use this 
parameter in a general classification of self-similar spherically symmetric perfect fluids. 
We note that E — Q for the flat FLRW solution, whereas the static solution has E <Q. 

Self-similar density perturbations of the flat FLRW solution can be studied 
covariantly following Ellis and B runi ||2^ . Th e fractional density gradient is 
characterized by a quantity L (see Appendix A.3). For the flat FLRW solution, L 
is zero, as expected. Orbits between the flat FLRW orbit and the V submanifold have 
L > and correspond to under-dense solutions, while orbits on the other side of the 
flat FLRW orbit have i < 0, corresponding to overdense solutions. 



5. Discussion 



In this paper we have studied timelike self-similar spherically symmetric perfect-fluid 
models using a dynamical systems approach. It continues the work initiated in a 
previous paper, dealing with the spatially self-similar models. By suitably adapting the 
variables to the homothetic Killing vector, the reduced phase space for these models 
becomes compact and regular, a desirable property that has never been obtained 
before. An advantage of a compact and regular phase space is that it enables a more 
thorough global picture of the space of solutions. Also, the current local treatment 
of the constraint makes it possible to avoid 'crushing' parts of the phase space. The 
solution space of the TSS models has a much more complicated structure than in the 
SSS case. In particular, the existence of a sonic line in one of the sonic surfaces has 
profound implications for the global structure. The sonic surfaces and the sonic line 
act as 'filters' determining differentiability of various solutions. 

To obtain a full global picture of self-similar spherically symmetric perfect-fluid 
models, it is necessary to match TSS and SSS regions. This has only been partially 
discussed here. A more complete discussion will be given in a subsequent paper. This 
will enable us to obtain a more detailed understanding of physical phenomena, such 
as formation of black holes and naked singularities. 

By allowing an additional coordinate dependence in the metric coefficients, 
the approach of this paper can easily be generalized to non-sclf-similar spherically 
symmetric models. This is of interest when studying critical behaviour in spherically 
symmetric gravitational collapse. However, the field equations will now be partial 
differential equations, and the phase space will no longer be compact. Nevertheless, 
for small perturbations, this generalization of the present approach may be useful. 
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Appendix A. Interesting physical quantities 

Appendix A.l. Fluid properties and Petrov type conditions 

The fluid expansion 9, the fluid shear a, the fluid acceleration scalar a and the Weyl 
scalar C are given by 

{u' + (1 - u^) [A +{2 + t+)u] } 9, 



V3(l-u2)3/2 

" " 3(1-^2)3/2 r' 

3V3(1-m2)3/2 



- {1-u^)(1 + 2J:+)u\9, 
\uu' + 3(1 - u'^)(uA ~ E+) I e, 



(A.l) 
(A.2) 
(A.3) 



3^/3 



|2SV -{l + q + S+)(l + 2E+) + KM^) 



where q was defined in (^^. Note that the magnetic part of the Weyl tensor is 
identically zero for all models. The spacetime is of Petrov type D if C 7^ 0, and of 
type if C = 0. 

Appendix A.2. Mass and energy 

For a general spherically symmetric spacetime the total mass-energy m between the 
centre distribution and some 2-space of symmetry is defined as (see, e.g., Misner and 
Sharp (29|, Hernandez and Misner ^ and Cahill and McVittie ||l|) 

V 2 



dR dR _K - {l + t^ 



A' 



R * " dx'^dx^ K ' ^^'^'^ 

where R — e^ D2{x) is the invariant radius of the symmetry surfaces. By definition, 
the matter is outside the gravitational radius 2m whenever ^ < 1 and inside the 
gravitational radius when ^ > 1. The latter case is associated with the existence of 
a black-hole apparent horizon or a cosmological apparent horizon at i? = 2m. 

It is possible to define the asymptotic energy per unit mass i? of a spherical shell 
as M 



1 




1 








2 


R 


^ 2 



(1 



Au) 



- 1 



where 



U 



A + u{l + T.+) 



kVT' 



(A.6) 



(A.7) 
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Appendix A. 3. Density perturbations 

A covariant approach to density perturbations has been given in [2^ . The fractional 
density gradient is defined as 

A, = = ft-^hydbfL. (A.8) 

In the reduced phase space variables, this leads to 

A, =0"e-*— ^(-u,l,0,0), (A.9) 

where 

^_ 2Au[^^ + (7-l)] +7[2»^ + (l + ^^^)S+] 



Appendix B. Coordinate and variable transformations 

Here we will give transformations to other coordinates and variables which have been 
used in the literature to study self-similar spherically symmetric models. The present 
variables lead to a compact and regular description everywhere in the TSS sector. 
This is not the case with previously used variables, as is easily seen in the equations 
below. The simple algebraic relations make it easy to identify the points, or even 
manifolds, where breakdowns occur. 



Appendix B.l. The variables used in the SSS region 

In a previous paper the form of the /i„ equation led us to compactify the 

variables with respect ioY — ^ 8'^ + 3-B| rather than 0. The resulting set of variables 
(Qoj Q+i Ci, v) are related to the variables used in this paper as follows: 

^ sgn(^) - _E+sgn(0") - isgn(^) _ _^ 



VTTk VTTk VTTk 

Note that since is a dominant quantity in the TSS region, Y is even more dominant. 
Thus we could have compactified the variables in the TSS region in the same way as in 
the SSS region with the additional replacement u — 1/v. However, the present choice 
of compactification result in simpler equations. The TSS variables are of course not 
compact in the SSS region, as 6 is not a dominant quantity there. This is easily seen 
in the above relations. 



Appendix B.2. Bogoyavlensky's variables 

Bogoyavlensky uses the homothetic approach, solves the constraint globally, and 
introduces the variables Q,w and u (as a starting point) Q. This results in an 
undesirable 'crushing' of the phase space. Bogoyavlensky's variables are related to 
the present ones by 

A 2T,^ , 
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Appendix B.3. The comoving (fluid) approach 

As the fluid 4-velocity singles out a preferred timelike direction in the spacetime, it 
should be beneficial to adapt the coordinates to this direction. The line element can 
be written as 

= -e^WdT^ + e^^^Ux^ + Y{\fT^dn^, (B.3) 

where A = T/X and dn'^ = d9^ + sin^ Odip'^. By defining T = exp[i - F{x)], 
X = exp[t — F{x) + a;], where the function F{x) satisfies 

dF ^ e^+^" 

the line element equation ( p3.3| ) can be written as 

ds^ = e^* {-Dldf + N^dx^ + Dldn^) (B.5) 

with 

i?? = (e* - , = "^1^^1+1^ = e-'^Y^. (B.6) 

Using that the fluid velocity is given by = 2^^*+^^, the above formulas can be 
written as 

" o*-2F /-I „,2\ f,j2 



with D2 the same as previously. The comoving approach has been used by many 
authors. The most useful variables produced so far to study self-similar spherically 
symmetric collapse of a perfect fluid in this approach are probably the ones obtained by 
Foglizzo and Henriksen |l^ (see also Bicknell and Henriksen ||^, |l^). Their variables 
{Np-H, fJ'F-HjVp-H} are related to the present ones by 

37S+ 371.(a--(1 + S+)' + A=) 

"-■'-id- \.(i-u')E,A ■ 

and vf-h = u. Here, we eliminated /it using the cons train t equation (H). Note the 



close connection between ^f-h and the mass function ( A. 5 ). Again it is obvious that 
these variables are not compact. However, they are useful for connecting the SSS and 
TSS regions and they are also very well suited for studying the phase space in the 
vicinity of the sonic line. 

Appendix B.^. The Schwarzschild approach 

One often represents a spherically symmetric line-element with its 'Schwarzschild' 
form: 

ds'^ ^~ FdT^ + GdR^ + R^dn^, (B . 9) 

where, in the self-similar case, F and G are functions of R/T only. By introducing 
R = e^B2^ , T = e*'e'^, where the function (j) satisfies the differential equation 

dx + 
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the line element ( |B.9| ) transforms into the diagonal homothetic form. The metric 
functions F and G become 

3il + t+)^e-^^0-^ ^ K 



F = 



G= : 2 —■ 



Defining = 1 — 2m/R recovers the definition of the mass function (A. 5). The 
radial 3- velocity of matter Vfi, i.e. the speed of the fluid with respect to the surfaces 
R — constant, is given by 

1 + 2j+ + uA 

The Schwarzschild approach has been used by several authors. The variables of Ori 
and Piran |^ are given by 

2 a jKu 
^ A + u{l + f^+) 

U = =r^. 

The orbits in the RC submanifold are parametrized by 

Do = lim —4^- (B.14) 

The sonic line is subsequently parametrized by the variable Y = M. /RD which on the 
sonic line is just 

. (B.15) 

2(7-l)+7S+ 

This approach was also used by Maison for studying the more general problem of 
non-universality in non-self-similar gravitational collapse. He used a set of variables 
(Am, Bm, Pm, vm)- In their renormalization group approach to the critical behaviour 
of spherically symmetric collapse, Hara et al [|o| used a set of variables {Nh, Ah,ojh) 
very similar to the variables of Maison (however, they solve the constraint globally). 
In the self-similar case, these sets of variables are given by 

^i, = A^i = li±^^— ^, Bl^N^'^ ^[ (B.16) 

UJH (1 - u') t+A 

Pm ^ -r- = i> ' = Vff = vr. 

Ah ^Ku 
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